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CUDA  Implementation  of  TE  -FDTD  Solution 
of  Maxwell’ s  Equations  in  Dispersive  Media 

Mohammad  Reza  Zunoubi,  Senior  Member,  IEEE ,  Jason  Payne,  and  William  P.  Roach 


Abstract — This  letter  presents  the  graphic  processor  unit  (GPU) 
implementation  of  the  finite-difference  time  domain  (FDTD) 
method  for  the  solution  of  the  two-dimensional  electromagnetic 
fields  inside  dispersive  media.  The  FDTD  is  truncated  by  the 
convolutional  perfectly  matched  layer  (CPML)  and  the  piece- 
wise-linear  recursive-convolution  (PLRC)  formulation  is  used  for 
modeling  dispersive  media.  By  using  the  newly  introduced  Com¬ 
pute  Unified  Device  Architecture  (CUDA)  technology,  we  illustrate 
the  efficacy  of  GPUs  in  accelerating  the  FDTD  computations  by 
achieving  significant  speedup  factors  with  great  ease  and  at  no 
extra  hardware/software  cost.  We  validate  our  approach  by  com¬ 
parison  to  exact  and  other  simulated  results,  which  show  favorable 
agreements.  The  effect  of  the  GPU-CPU  memory  transfers  on  the 
speedup  factor  will  be  also  studied. 

Index  Terms — Convolutional  perfectly  matched  layer  (CPML), 
dispersive,  finite-difference  time  domain  (FDTD),  graphic  pro¬ 
cessor  unit  (GPU). 


I.  Introduction 

THE  interaction  of  electromagnetic  fields  with  human  tis¬ 
sues  and  the  resulting  biological  effects  have  been  of  great 
interest  to  both  commercial  and  government  sectors.  The  finite- 
difference  time  domain  (FDTD)  technique  has  been  a  method 
of  choice  for  many  years  in  studying  this  interaction.  FDTD’s 
popularity  is  partly  due  to  its  straightforward  implementation 
and  its  ability  to  model  a  broad  range  of  exposure  conditions 
for  anisotropic,  nonlinear,  and  dispersive  media.  However,  the 
FDTD  algorithm  can  be  computationally  expensive,  and  signif¬ 
icant  run-times  may  exist  for  certain  applications. 

Recently,  NVIDIA  has  introduced  a  general-purpose  parallel 
computing  architecture  called  Compute  Unified  Device  Archi¬ 
tecture  (CUDA)  that  leverages  the  parallel  compute  engine  in 
NVIDIA  graphic  processor  units  (GPUs)  to  solve  many  com¬ 
plex  computational  problems  in  a  fraction  of  the  time  required 
on  a  CPU  [1].  Through  CUDA,  the  vast  variety  of  GPU  memory 
hierarchy  may  be  accessed,  which  allows  for  maximum  opti¬ 
mization  of  computational  schemes  that  are  parallel  in  nature, 
such  as  FDTD. 

As  commercial  sectors  like  Remcom,  Computer  Simulation 
Technology  (CST),  Acceleware,  Schmid  &  Partner  Engineering 
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AG  (SPEAG),  and  Agilent  have  integrated  the  GPU-acceler¬ 
ated  FDTD  technique  into  their  software  packages  to  treat  elec¬ 
tromagnetic  problems  in  general  dispersive  and  nondispersive 
media,  researchers  in  academia  have  also  focused  on  taking  ad¬ 
vantage  of  inherent  attributes  of  the  GPUs  for  parallel  com¬ 
puting.  Such  implementations  that  use  the  CUDA  technology 
include  the  works  presented  in  [2]  and  [3]  originally  and  the 
research  presented  by  Sypek  et  al.  for  a  TMZ  solution  of  elec¬ 
tromagnetic  fields  [4],  the  double-precision  implementation  of 
the  FDTD  on  the  Tesla  GPU  [5],  and  the  work  reported  in  [6] 
more  recently,  although  earlier  successful  FDTD  GPU-acceler¬ 
ated  methodologies  using  less  versatile  APIs  such  as  Brook  or 
OpenGL  were  reported  starting  a  few  years  ago. 

Here,  we  present  the  implementation  of  the  two-dimensional 
(2-D)  FDTD  method  in  dispersive  media  on  GPUs  using  CUDA 
technology.  The  convolutional  perfectly  matched  layer  (CPML) 
is  used  to  truncate  the  solution  domain  [7],  and  a  plane-wave  in¬ 
cident  field  is  utilized  to  simulate  realistic  exposure  scenarios. 
The  software  validity  is  established  by  comparison  to  exact  re¬ 
sults  for  the  reflection  coefficient  of  a  dispersive  half-space  and 
other  simulated  results  for  more  complex  geometries.  Speedup 
factors  are  reported  that  show  promise  for  the  rapid  solution  of 
large-scale  electromagnetic  radiation  and  scattering  problems. 
We  also  study  the  effects  of  the  memory  transfer  between  the 
GPU  and  CPU  on  the  speedup  factor. 

II.  FDTD  Formulation 

In  the  following,  we  present  the  analysis  of  TEZ -polarized 
plane  wave  penetration  through  two-dimensional  dispersive 
media.  As  such,  we  follow  the  procedure  outlined  in  [7]  for 
the  piecewise-linear  recursive-convolution  (PLRC)  method  for 
the  single-pole  Debye  representation  of  the  dispersive  media 
along  with  the  CPML  boundary  condition,  which  is  not  only 
efficient  in  minimizing  the  memory  requirements,  but  also 
the  most  accurate  form  of  absorbing  material.  This,  in  turn, 
yields  the  following  form  of  the  FDTD  update  equation  for  the 
rr-component  of  the  electric  field: 

Tp  in+1  _  „  |  Tp  in  i  | 

•r  U+ 1  / 2, y  —  CaU+l/2,i^|  j+l/2,y  cb\i+l/2,j 

Tf  In  _  TT  In 

Jp^\i+l/2J+l/2  JppU+l/pj-l/p 
Kyj 

+  ^EX,y\i+1/2,j) 

+  Cc\i+l/2,j^Ex\’t+1/2j 

where  Ca ,  C*,,  and  Cc  are  the  coefficients  associated  with 
space-time  discretization  and  material  properties,  nyj  is  the 
scaled  tensor  parameter,  i>Ex,y  is  the  auxiliary  array  for  the 
CPML,  and  v/jex  is  the  recursive  accumulator  variable  given 
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in  [8]  and  defined  in  terms  of  eS9  ,  and  r,  which  are  the 
static  relative  permittivity,  the  relative  permittivity  at  infinite 
frequency,  and  relaxation  time,  respectively,  for  the  Debye 
material. 

To  model  the  tissue  material,  we  generate  an  id  file  that  is 
used  to  assign  the  constitutive  parameters  to  the  corresponding 
field  components.  For  normal  materials,  and  er,  pr  and  a  are 
assigned  and  used  to  calculate  the  regular  FDTD  updating  coef¬ 
ficients.  For  dispersive  materials,  the  Debye  parameters  £oo,  eS9 
and  r  are  assigned  and  used  to  calculate  the  dispersive  FDTD 
updating  coefficients  for  the  electric  fields  and  for  the  recursive 
accumulator  variable  . 

III.  CUDA  IMPLEMENTATION 

A.  Memory  Management 

Since  we  are  implementing  a  TEZ  formulation,  we  allocate 
and  initialize  on  the  global  memory  of  the  device  (GPU)  three 
one-dimensional  arrays  for  the  Ex,  Ey,  and  Hz  components  of 
the  fields,  four  auxiliary  ^-arrays  for  the  CPML  such  as  ^  , 

two  arrays  for  recursive  accumulator  variables  and  ipEy , 
and  two  arrays  for  the  incident  electric  and  magnetic  field  com¬ 
ponents.  Note  that  all  these  arrays  will  have  to  be  updated  during 
the  FDTD  time-marching  process.  However,  once  on  the  device , 
we  will  use  the  fast  but  limited  shared  memory  as  described  in 
the  next  section  to  handle  updating  equations.  Additionally,  we 
need  to  allocate  and  assign  two  integer  arrays  tissuelDx  and  tis- 
suelDy  that  contain  the  indices  assigned  to  the  different  tissue 
material  under  exposure,  and  16  arrays  for  coefficients  ca, 
and  cc,  scaled  tensor  parameters  ft’s,  and  the  constants  used  in 
calculating  auxiliary  ^-arrays.  However,  since  no  updating  of 
these  arrays  is  needed  during  the  time-marching  process,  we  use 
the  texture  memory  of  the  device  for  storage,  which  allows  for 
fast  read-access  when  updating  the  field  quantities. 

B.  Device  Kernels 

For  efficient  implementation  of  our  FDTD  approach,  we  di¬ 
vided  the  device  implementation  into  two  kernels :  one  to  update 
the  magnetic  field  component  and  another  to  update  the  electric 
field  intensities.  To  take  full  advantage  of  the  single-instruction 
multiple-date  (SIMD)  architecture  of  the  GPU  and  overcome 
the  memory  latency,  we  group  the  threads  into  i  x  j  2-D  thread 
blocks.  One  thread  is  assigned  to  an  output  element  while  al¬ 
lowing  four  halo  regions  to  include  the  neighboring  elements 
that  are  needed  in  the  FDTD  updating  equations.  This  method 
allows  the  threads  within  each  small  block  to  access  their  own 
shared  memory,  the  latency  of  which  is  two  orders  of  magni¬ 
tude  lower  than  that  of  global  memory.  Therefore,  during  each 
time  step,  the  field  and  ^-arrays  are  loaded  into  the  shared 
memory,  the  updating  equations  are  operated  upon,  and  the  re¬ 
sulting  updated  electric  and  magnetic  fields  are  stored  in  the 
global  memory  to  be  read  by  the  host. 

Note  that  the  thread  dimensions  i  and  j  of  each  block  need  to 
be  optimized  in  order  to  take  full  advantage  of  the  GPU  compu¬ 
tational  capabilities.  This  will  be  addressed  in  the  Section  IV. 

Following,  we  illustrate  our  CUDA  implementation  by 
describing  the  kernels  on  the  device  side  that  use  the  texture , 
shared ,  and  global  memory  hierarchies  to  perform  the  FDTD 
updating  equations  efficiently.  Note  that  before  invoking  these 
kernels,  the  main  program  on  the  host  side  should  manage  the 


device  memory  allocations,  copying  data  from  host  to  device 
and  copying  the  data  after  FDTD  updating  (invocation  of 
kernels)  from  the  device  to  the  host : 

•  Load  the  read-only  data  in  the  texture  memory: 

texture<int,  1,  cudaReadModeElementType>ttissueIDx; 

•  Repeat  for  all  other  constant  arrays; 

•  update  H<Grid  of  thread  blocks,  Size  of  thread  block> 

{ 

Compute  spatial  coordinate  of  current  thread 
Load  Ex  and  Ey  from  the  global  memory  into  the 
shared  memory  to  avoid  GPU’s  memory  latency; 
Perform  FDTD  updating  equations  for  H  while 
reading  constant  arrays  from  the  fast  texture  memory; 
Copy  data  from  the  device  global  memory  back  to 
the  host’, 

} 

•  update  E<Grid  of  thread  blocks,  Size  of  thread  block> 

{ 

Compute  spatial  coordinate  of  current  thread 
Load  Hz  from  the  global  memory  into  the  shared 
memory  to  avoid  GPU’s  memory  latency; 

Perform  FDTD  updating  equations  for  E  while 
reading  constant  arrays  from  the  fast  texture  memory: 
—  Compute  CPML  modifiers  (d_PsiExy,  d_PsiEyx) 
—  Perform  updating  of  E-fields  including  the  CPML 
for  no-dispersive  regions 
—  Update  ipEx  and  ^ev  matrices  for  dispersive 
material 

—  Update  the  E  fields  in  dispersive  regions 
—  Correct  for  plane  wave  interface 
Copy  data  from  the  device  global  memory  back  to 
the  host ; 

} 

C.  Arithmetic  Instructions 

To  maximize  the  computational  efficiency  of  our  kernels ,  we 
use  the mul24  and / divide  f  functions  for  integer  multi¬ 
plications  and  floating-point  divisions,  respectively. _ mul24 

provides  24-bit  integer  multiplication  in  four  clock  cycles 
compared  to  the  16  clock  cycles  for  the  conventional  32-bit 

multiplication,  and  / divide/  performs  single-precision 

floating-point  division  in  20  clock  cycles,  which  is  superior  to 
the  36  clock  cycles  typically  needed  for  dividing  floating  point 
values  [1]. 

IV.  Results 

Before  studying  the  advantages  of  the  CUDA  implementa¬ 
tion  of  the  FDTD  for  dispersive  media,  we  demonstrate  the  ac¬ 
curacy  of  our  developed  tool  by  first  analyzing  the  problem  of  a 
Gaussian  plane  wave  propagating  in  vacuum  normally  incident 
on  a  vacuum-water  interface.  The  problem  space  is  divided  into 
1024  cells  in  both  x-  and  ^-directions  with  a  cell  size  of  37.5  pm, 
which  ensures  the  accuracy  of  the  results  up  to  400  GHz  in 
free-space  regions  based  on  the  20-cells/A  criterion.  The  time 
step  is  chosen  according  to  the  Courant  limit  as  At  =  0.0884  ps, 
and  we  use  an  eight-cell  CPML  absorbing  boundary  to  trun¬ 
cate  the  volume.  To  model  the  frequency  dependence  of  water, 
the  Debye  parameters  of  £s  =  81,  £oo  =  1.8,  a s  =  0,  and 
r  =  9.4e  —  12  are  used.  The  FDTD  calculation  is  performed 
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Fig.  1.  The  magnitude  of  the  reflection  coefficient  calculated  at  the  vacuum- 
water  interface. 


y 


(a)  (b) 

Fig.  2.  (a)  Problem  geometry  for  a  combination  of  dispersive  and  nondis- 
persive  materials,  (b)  Magnitude  of  the  reflection  coefficient  calculated  at  the 
vacuum-Debye  interface. 


for  1300  time  steps,  and  results  are  recorded  at  the  cell  loca¬ 
tion  (504,  504),  which  is  on  the  vacuum-water  boundary.  The 
incident  field  due  to  the  Gaussian  pulse,  which  has  duration  of 
256  time  steps,  is  subtracted  from  the  total  field  Ey  to  obtain  the 
reflected  field  due  to  the  material  discontinuity.  The  reflection 
coefficient  versus  frequency  is  then  calculated  by  normalizing 
the  Fourier  transform  of  the  reflected  field  by  the  Fourier  trans¬ 
form  of  the  incident  pulse.  Results  can  be  seen  in  Fig.  1,  where 
a  comparison  is  made  with  the  exact  solution,  illustrating  the 
adequacy  and  accuracy  of  our  FDTD-CUDA  implementation. 

To  further  demonstrate  the  accuracy  of  our  implementation, 
we  consider  a  2-D  region  shown  in  Fig.  2(a),  which  consists 
of  layers  of  vacuum,  a  Debye  slab,  vacuum,  a  dielectric  slab, 
and  vacuum  again,  and  we  compute  the  reflection  coefficient  as 
discussed  on  the  vacuum-Debye  material  interface  versus  fre¬ 
quency.  Note  that  all  the  dimensions  are  in  terms  of  the  FDTD 
voxels.  The  Debye  material  properties  are  for  water  as  indi¬ 
cated  previously,  and  the  dielectric  slab  has  sr  =  2.  Fig.  2(b) 
shows  the  results  along  with  comparisons  to  the  corresponding 
results  obtained  from  the  Remcom  XFDTD  software,  where  a 
very  good  agreement  is  observed.  Note  that  for  the  XFDTD  re¬ 
sults,  the  region  is  extended  in  the  ^-direction  by  400  FDTD 
cells  to  mimic  a  2-D  problem.  In  our  computation,  the  same 
spatial  discretization  as  the  above  example  was  chosen  to  assure 
the  accuracy  of  the  results  up  to  400  GHz  while  the  simulation 
was  performed  for  1100  time  steps,  which  insured  achieving  the 
steady-state  condition. 


18000  - 

16000  -  / 


Time  (ns) 

(a)  (b) 

Fig.  3.  (a)  Cuvette  filled  with  water  modeled  as  Debye  material,  (b)  Incident 
and  total  electric  field  calculated  at  the  center  of  the  cuvette. 

We  next  demonstrate  the  applicability  of  our  tool  to  cal¬ 
culating  the  absorbed  dose  of  dispersive  media  exposed  to 
nonionizing  ultrawideband  (UWB)  electromagnetic  pulses  of 
nanosecond  duration.  This  dosimetry  is  critical  for  establishing 
dose-response  curves  for  nanosecond  electromagnetic  pulses. 
Here,  we  consider  the  problem  of  exposing  water  inside  a 
cuvette  to  a  short  pulse  described  as  [9] 

Ey  =  Eo(e~at  -  e~'dt)  with  Eq  =  18.5  KVm-1, 

a  =  108  s-1,  0  =  2.0  X  1010  s_1  (2) 

which  has  a  rise  time  of  150  ps  and  a  width  of  10  ns.  The  cu¬ 
vette  is  1  x  1  cm2  with  1-mm-thick  plastic  walls,  as  shown  in 
Fig.  3(a),  and  is  filled  with  water  that  has  the  dispersive  proper¬ 
ties  described.  Since  the  penetration  of  the  electric  field  compo¬ 
nent  in  the  direction  of  polarization  is  defined  by  the  rise  time 
and  pulsewidth,  we  plot  both  the  incident  field  and  the  total  field 
inside  water  in  the  middle  of  the  cuvette  in  Fig.  3(b).  This  plot 
indicates  that  the  pulse  inside  water  is  a  superposition  of  a  short 
pulse,  induced  during  a  fast  rise  time,  and  recursively  longer 
pulses  induced  by  the  slow  variation  in  the  incident  pulse.  This 
is  similar  to  the  observation  made  in  [9] . 

As  discussed  in  the  previous  section,  the  size  of  the  thread 
blocks  loaded  into  shared  memory  could  have  a  significant  ef¬ 
fect  on  the  performance  of  the  GPUs.  Accordingly,  a  study  was 
performed  to  identify  the  optimum  block  dimension  for  our 
FDTD  computations.  We  executed  the  CUDA  program  for  2000 
time  steps  for  various  thread  arrangements  following  the  format 
of  ( ixj )  =  (2xj),  (4 x  j),  ( 8x7 ),  (16 x  j),  (32 x  j),  (64 x  j), 
(128  x  j),  and  (256  x  j )  while  setting  j  =  2,  4,  8,  16,  32,  64, 
128,  and  256  with  a  problem  space  of  size  1024  x  1024  FDTD 
cells.  This  study  indicated  that  a  block  size  of  ixj  =  16  x  4, 
which  results  in  64  threads  per  block  and  uses  256  bytes  of  the 
available  16  kB  of  memory,  yields  the  optimum  performance 
with  a  GPU  computation  time  of  7.87  s. 

Having  identified  the  optimum  block  dimension,  we  next  in¬ 
vestigate  the  speedup  factors  achieved  using  the  GPU  as  we  in¬ 
crease  our  FDTD  domain  size  from  128  x  128  cells  (16  384  grid 
points)  to  2048  x  2048  cells  (4  194  304  grid  points).  The  FDTD 
computations  are  performed  for  2000  time  steps  on  both  C  and 
CUDA  versions  of  our  code.  The  C  implementation  was  per¬ 
formed  on  Intel  Core2  Quad  Q6600,  2.40  GHz,  and  3.25  GB 
of  RAM.  The  CUDA  tool  was  run  on  an  NVIDIA  GeForce 
9800  GT,  which  has  112  stream  processors,  1024  MB  standard 
memory,  a  memory  bandwidth  of  57.6  GB/s,  and  a  shader  pro¬ 
cessing  rate  of  504  gigaflops.  The  computation  times  for  both 


ZUNOUBI  et  al.:  CUDA  IMPLEMENTATION  OF  TE2-FDTD  SOLUTION  OF  MAXWELL’S  EQUATIONS  IN  DISPERSIVE  MEDIA 


759 


TABLE  I 

CPU  and  GPU  FDTD  Computation  Time 


Grid  points 

CPU  (s) 

GPU  (s) 

Speedup 

16,384 

9.28 

0.39 

23.8 

65,536 

43.81 

0.76 

57.6 

262,14  4 

186.4 

2.20 

84.72 

1,048,576 

794.5 

8.05 

98.70 

4,194,304 

3114.7 

31.6 

98.56 

C  and  CUDA  implementations  as  well  as  the  achieved  speedup 
factors  are  given  in  Table  I.  These  results  show  a  significant  im¬ 
provement  in  the  FDTD  calculations  time  for  the  GPU  versus 
the  CPU  implementation,  and  it  is  seen  that  speedup  factor  of 
98.56  can  be  obtained  for  up  to  4  194  304  grid  points.  It  is  ex¬ 
pected,  however,  that  with  the  newly  manufactured  NVIDIA 
Tesla  Cl 060  GPUs,  these  computations  can  be  carried  out  for 
hundreds  of  millions  of  FDTD  grid  points  with  even  greater 
computational  efficiency.  It  is  also  to  be  noted  that  the  time 
taken  to  transfer  the  data  from  the  host  to  the  device  and  vice 
versa  was  negligible  and  was  found  to  be  only  0.17  s  for  the 
4  194  304  FDTD  cells  when  transferring  the  initial  values  of  the 
fields  and  material  arrays  to  the  device. 

It  is  known  that  for  some  practical  applications,  there  is  a  need 
for  recording  the  time-domain  field  values  at  several  locations 
inside  the  FDTD  solution  domain  during  the  time-marching 
process.  When  working  with  the  GPU  kernels ,  there  are  three 
possible  approaches.  One  is  to  store  such  data  on  the  device  as 
one  iterates  through  all  time  steps  and  transfer  the  data  to  the  host 
at  the  end  of  the  time-marching  process.  This  would  obviously 
have  less  effect  on  the  speedup  factors  that  can  be  achieved  by 
performing  the  computations  on  the  GPU.  However,  this  would 
require  allocating  and  using  part  of  the  global  memory  of  the 
device ,  further  limiting  the  number  of  unknowns  that  can  be 
treated  using  one  single  GPU.  The  second  option  would  be  to 
save  the  sample  points  on  the  device  only  during  each  single  time 
step  and  transfer  them  to  the  host  after  the  updating  equations  are 
performed  to  then  be  saved  in  their  corresponding  arrays  on  the 
host.  This  approach  needs  a  small  amount  of  the  device  memory 
usage  and  should  not  affect  the  speedup  factors  considerably.  The 
third  method  would  be  to  transfer  the  entire  array  of  the  desired 
field  component(s)  to  the  host  after  every  time  step,  and  then  save 
the  required  sample  point  locations  in  their  corresponding  arrays 
on  the  host ,  thus  eliminating  any  need  for  additional  storage  on  the 
device.  By  eliminating  the  memory  requirements  on  the  device , 
all  the  device  memory  space  can  be  devoted  to  performing  the 
FDTD  updating  equations  at  the  cost  of  lowering  speedup  factors 
due  to  the  memory  transfer  at  each  time  step  from  the  device  to 
the  host  specially  as  the  problem  size  grows. 

In  order  to  investigate  the  data-transfer  effect,  we  studied 
the  effect  of  sampling  the  ^/-component  of  the  electric  field  at 
10  equally  spaced  locations  along  the  x-axis  at  the  middle  of 
the  FDTD  solution  domain  for  8192  time  steps,  first  by  com¬ 
puting/saving  the  sampled  field  values  on  the  device  for  each  time 
step  and  then  transferring  the  data  onto  the  host ,  and  second,  by 
transferring  the  whole  array  of  the  ^/-component  of  the  electric 
field  after  each  time  step  and  storing  only  the  sample  points  results 
on  the  host.  As  expected,  for  the  first  approach,  no  appreciable 


reduction  in  the  speedup  factors  was  observed,  while  some  device 
memory  needed  to  be  allocated  for  the  storage  of  the  sample 
points.  Following  the  second  method  reduced  the  speedup  factors 
from  98.7  and  98.56  s  for  1  048  576  and  4  194  304  unknowns  to 
39.5  and  43.5  s,  respectively,  further  demonstrating  the  fact  that 
transferring  large  arrays  of  field  components  from  the  device  to 
the  host  can  reduce  the  speedup  factors  considerably. 

V.  Conclusion 

We  presented  a  CUDA-based  FDTD  analysis  of  dispersive 
media  on  an  NVIDIA  9800  GT  GPU  card.  The  accuracy  of  our 
implementation  was  illustrated  by  comparing  our  results  with 
their  corresponding  analytical  and  simulated  results.  Very  good 
agreement  was  observed.  We  further  investigated  the  computa¬ 
tional  advantage  of  performing  the  FDTD  equations  on  GPUs 
by  solving  for  the  electromagnetic  fields  in  dispersive  media 
for  a  problem  space  of  2048  x  2048  voxels.  It  was  shown  that 
speedup  factors  of  up  to  98.56  can  be  achieved  by  facilitating  the 
shared  memory  of  the  thread  blocks  and  by  using  the  optimized 
block  dimensions  and  arithmetic  of  the  CUDA  technology.  The 
effects  of  memory  transfer  between  the  device  and  the  host  were 
also  studied  for  10  sample  locations,  illustrating  the  fact  that  if 
large  arrays  of  data  are  to  be  transferred  from  the  device  to  the 
host  at  each  time  step,  considerable  reduction  in  speedup  factors 
may  occur.  We  are  currently  implementing  the  more  practical 
three-dimensional  (3-D)  version  of  our  tool  and  expect  to  ob¬ 
tain  considerable  speedup  factors  with  minimal  or  no  additional 
software/hardware  resources,  although  it  is  known  that  the  3-D 
implementation  would  not  yield  large  speedup  factors  achieved 
with  the  2-D  tool  since  the  GPU  structure  is  optimally  made  for 
the  support  of  the  2-D  images. 
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